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Abstract 

The development of reliable methods for estimating microcanonical averages constitutes an im- 
portant issue in statistical mechanics. One possibility consists of calculating a given microcanonical 
quantity by means of typical relations in the grand-canonical ensemble. But given that distinct 
ensembles are equivalent only at the thermodynamic limit, a natural question is if finite size ef- 
fects would prevent such procedure. In this work we investigate thoroughly this query in different 
systems yielding first and second order phase transitions. Our study is carried out from the direct 
comparison with the thermodynamic relation (§§), where the entropy is obtained from the den- 
sity of states. A systematic analysis for finite sizes is undertaken. We find that, although results 
become inequivalent for extreme low system sizes, the equivalence holds true for rather small L's. 
Therefore direct, simple (when compared with other well established approaches) and very precise 
microcanonical quantities can be obtained from the proposed method. 

PACS numbers: 05.10.Ln, 05.20.Gg, 05.50.+q 
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I. INTRODUCTION 



The development of efficient Monte Carlo (MC) methods constitutes a key problem in 
statistical mechanics. Typically, numerical simulations are performed by one-flip algorithms 
that generate a grand-canonical ensemble, in which intensive quantities are held fixed 
However, in some cases, a different ensemble would be more appropriate. For example, 
strong first-order phase transitions become extremely hard to simulate by one-flip grand- 
canonical schemes: the presence of different phases separated by large free-energy barriers 
makes the system to be trapped into metastable states at the phase coexistence, even for 
small system sizes. Although different grand-canonical procedures have been proposed 2^ 
4|, the microcanonical ensemble is also an appropriate way to circumvent these problems. 
In this case, the simulation is carried out for fixed energies and intensive thermodynamic 
quantities are treated as external variables, hence avoiding large entropic barriers. 

Different microcanonical schemes have been proposed in the last years. Entropic sampling 
5], broad histogram method js] and the Wang-Landau (WL) method [3] are some examples 
of procedures in which the density of states (DOS) is estimated numerically. Other micro- 
canonical approaches, not requiring the knowledge of the DOS, have also been developed. 
In such cases, the temperature and other intensive quantities are obtained from auxiliary 
relations. For example, Creutz j§| has generated a microcanonical ensemble by assuming a 
canonical distribution of the energies carried by a "demon" , where the total energy is not 
strictly conserved, but it fluctuates above a constant lower bound. More recently, Martin- 
Mayor j^l proposed a method where the temperature is obtained from an ensemble of fixed 
energy (including the kinetic energy). The temperature is calculated from the fluctuations 
of the spin part. 

To analyze the discontinuous transitions in the Blume- Emery-Griffiths (BEG) jic| and 



asymmetric Ising models [ll|, it has been employed a "trick" [12j, where intensive quanti- 
ties are calculated through expressions originally derived in the grand-canonical ensemble 
Although the equivalence is granted in the thermodynamic limit, Gibbs ensembles may be in 



equivalent for finite systems, including short range interactions at the phase coexistence [13]. 
Therefore, in order to the protocol can be extended to more general situations, one should 
test under what conditions the finite size effects hinder the equivalence between methods. 
In other words, it should be verified that intensive quantity as a matter of fact correspond 



to the genuine microcanonical temperature (obtained from the derivation of entropy with 
respect to the energy). 

The first goal of this paper is then to answer the above query. For so, different aspects 
of the method will be exemplified by means of distinct lattice systems yielding first and 
second-order phase transitions. We first address the Ising model, for which the DOS are 
known exactly fl^j ]. Then we consider as next examples the Blume-Capel (BC) and the 
Potts models. Although they do not present exact DOS, we are going to compare with the 
very efficient WL sampling as a benchmark. The BC model is an interesting case, since 
its DOS has been obtained by performing a random walk in the space of two parameters 



(in similarity to several lattice models presenting distinct particle interactions) 111, |l5] . We 



intend to verify if the calculation of the temperature will be changed by different restrictions 



in the random walk (as performed in Ref. 15|). 



The Potts models is also a very interesting test. Unlike the previous cases, its discontinu- 
ous transitions (yielding for q > 4) presents a genuine microcanonical feature, the existence 



of a loop 



3,y, 



16| . Thus, it is important to verify if our approach not only reproduces this 



remarkable signature but also is equivalent to the (|§). In addition, we are also exploit- 
ing the interesting q = 4 case, that although yielding a continuous transitions, it possesses 
distinct behaviors inclu ding logarithmic scaling corrections ji], 0, Il8| and a double peak 



probability distribution 19j. As it will be shown further, the methods becomes equivalent 
in all above models for relatively small system sizes L. The equivalence includes not only the 
temperatures but also all extrapolated thermodynamic limit points. However, for extreme 
small L's, the results do not agree and thus the intensive quantity can not be recognized 
as the thermodynamic temperature. We will present a detailed analysis showing how the 
methods converge when L increases. 

A second contribution here is to exploit the advantages. Besides its low computational 
cost, it does not require criterion for achieving convergence of results. Other immediate 
advantage is that intensive quantities are evaluated directly from standard numerical sim- 
ulations and become more precise as the system size increases. In addition, the method is 
very easily extended for other lattice models. 

This paper is organized as follows: In Sec. II we review the methods for calculating the 
temperature. In Sec. Ill we show the numerical results for the models and in Sec. IV we 
present our conclusions. 



II. MICROCANONICAL TEMPERATURE 



For a given system size L and energy per site e = E/V (where V = L d and d is the 
dimension), the inverse of microcanonical temperature 0l is obtained through the expression 

Pl = (1) 

where s = s(e,L) = |f lnfl(E, L) is the entropy per site and /3l = l//csT £ . The quantity 
Q(E, L) denotes the DOS for given E and L. For the Ising model, the Q(E, L) is known 
exactly [l^]. For the other models, we shall estimate Q(E,L) using the WL sampling 

The WL sampling is a powerful technique to calculate Q(E, L) by carrying out a random 
walk in energy space with an acceptance probability proportional to 1/Q(E, L), i.e., 

where Ei and Ej are the energies of the current and a possible new configuration, respectively. 
For each new accepted configuration an energy histogram H(E) is accumulated. 

During the random walk, whenever a move to a configuration with energy E is accepted, 
Q(E, L) is updated by multiplying it by a "modification factor" / > 1 that accelerates the 
diffusion of the random walk, and an unit is added to the histogram H(E). The initial choice 
of / is f = e = 2.71828 .... Q(E, L) is multiplied by / until the accumulated histogram 
H(E) becomes flat. We then reduce / by setting / — > y/J, and resetting H(E) = for all 
energy values. The simulation converges to the true value of Q(E, L) when / approximates 
to the unit. In particular, in this work we use the improved Wang-Landau sampling proposed 



by Cunha-Netto et al. [20[. Their approach use adaptive energy windows to eliminate border 
effects that affect the density of states mainly of q-states Potts model, which is our case here. 
In our simulations the criterion of flatness was taken as each value of the histogram reaching 
at least 80% of the mean value (H(E)) for the BC model and 90% for the Potts model. The 
histograms are generally checked after each 10000 MC steps. Here we performed 10 different 
runs for the same L with different initial seeds in order to reduce statistical fluctuations. 
For the BC model, the DOS is obtained by performing a random walk for two parameters 
E and E 2 (E 2 = ^ °f )• An immediate advantage of the WL is that a single run gives the 
DOS for the whole range of energy, which provides the calculation of canonical averages for 
any temperature. 



Now, following Ref. 12] we proceed to obtain the temperature T with respect to the 
microcanonical ensemble. The method consists in writing down the probabilities of differ- 
ent microscopic configurations in the grand-canonical ensemble and further resorting the 
equivalence of ensembles. 

The probability distribution P(er) of a microscopic configuration a = (<7i, 02, cry) in 
the grand-canonical ensemble is given by P(cr) = exp{— /3"H(cr)}/S, with Hamiltonian H{o) 
reading 

U(a) = -J^attTj, (3) 



(id) 



for the Ising model and 



for the BC model and 



W(<r) = -J5>*i + A XH (4) 



for the Potts model. Parameters J and A are the energy between two nearest-neighbor spins 
and the crystalline field, respectively. The spin variable a takes the values —1 or +1 for the 
Ising model, —1,0 or +1 for the Blume-Capel and 0, 1, q — 1, for the Potts model. In all 
cases, the summations are restricted over nearest neighbor sites. 

By considering the transition — 1 -<->• 1 for Ising and BC models and denoting a k by a 
microscopic configuration which differs from a only by the value of the spin at the site k, that 
is, a k = (o"i, <72, — .., cry), the ratio between P(er) and P(cr fc ) in the grand-canonical 
ensemble is given by 

-^ = exp{2M[0 fc (a)]}, (6) 

where 

(f)k(cr) = J^2a k+5 , (7) 

5 

whose summation is performed over 5 nearest neighbor sites. If the average of an arbitrary 
state function in the grand-canonical ensemble is given by ( f{o)) gc = Yla f( a )P( a )i f rom 
Eq. ()6]) we have that 

</(*)>* = (iV) exp{2(3a k lMv)}}} 9 c. (8) 
By taking Eq. (jSJ) for f(a) given by f(a) = 5(<Tk, +l)8(<f>k(a), E), we have that 

-2pE = (^fc,~lWfc(>),P))gc (q) 

(S(a k , +l)6{M<r), E)) 9 c 1 J 
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where E denotes one of all possible values of <fi k {(j). In the microcanonical ensemble, the 
energy per site e = E/V is held fixed. Assuming the equivalence between the grand-canonical 
and microcanonical ensembles, we get the following expression 

which allows us to obtain the temperature T with respect to the microcanonical ensemble. 
A similar procedure can be performed for the Potts model. By choosing two particular 
states k* and k** (ranging from to q — 1) with respective transition k* <r> k** and the state 
function h(a) = S(a k , k*)S((f) k (a), E) we have, by appealing to the equivalence of ensembles, 
that 

(11) 



-PE _ (S(ak,k*)5((f) k (a),E)} r 



(5(a k ,k**)5(M°),E)) mc ' 
where 4> k (cr) is given by (p k (cr) = JJ2s(^k**,a k+ s — 8k*,a k +s)- Since the above formulae does 
not specify the dynamics, they are valid for different classes of microcanonical algorithms. 



III. NUMERICAL RESULTS 

Numerical simulations have been performed in a square lattice with L 2 sites. The micro- 
canonical dynamics is composed of two parts. In the first part, a given site of the lattice is 
randomly chosen and its spin is changed to one of its all possible values. In the second part, 
two sites of the lattice, also randomly chosen, have theirs spins interchanged. The above 



dynamics are accepted only if the total energy remains unchanged. It is wort 



i mentioning 



that the actual MC algorithm is quite different from those studied in Refs. 12] , where both 
energy and magnetization are strictly conserved. Here the particle moves are accepted only 
when the total energy does not change. The number of species (spins) is not necessarily 
conserved. 

By applying the logarithm on both sides of Eq. ( flOl) . we get the following expression 
(written in units of J and ks) 

\KRn = ~Y, (12) 

where the right-hand of Eq. ffTU]) was written as R n . The possible values of the quantity 
E = 4> k (cr) = JJ2s a k+s are given by nJ, where n takes the values n = —4,-2,0,2,4, 
for the Ising model and n = —4, —3, 3,4, for the BC model. Thus, from the above, by 
calculating \nR n numerically for all possible values of 4> k (cr), the temperature is extracted 
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from the inverse of the slope of Eq. (fT2l) . A similar procedure is done for the Potts model, 
where we have 

71 

\nR n = --, (13) 

where n assumes the values —4, —3, 3, 4 for all values of q. 

In the first analysis, we study the validity of Eqs. (|T2l) and ( TLB"]) for different L, as showed 
in Fig. []Ja) and (b) for the Ising model. Continuous lines and symbols denote standard 
and present results, respectively and temperatures calculated for L = 4 are exact in both 
cases. Comparison between intensive quantities show that they are slightly different for the 
smallest L's (L = 4 and L = 6). Inspection of part (b) reveals us that in these cases, the 
quantity lni? n is not linear in n, and hence Eq. ffT2l) does not hold. The non validity of 
Eq. ( fT2l) for extreme small L is exemplified by evaluating the numerator and denominator 
of Eq. ffTUl) for e = —1.25 and L = 4. Since the number of configurations is small, both 
quantities are zero for n = 0, whereas for n = 4 (n = —4) the numerator (denominator) is 
null. By increasing L the number of configurations becomes large in such a way that the 
linear dependence between lni? n and n is achieved. Only in this regime we can evaluate T 
from Eq. ( JT2l) . In practice, estimates become equivalent for rather small system sizes. For 
example, for the Ising model the difference between estimates is in the third decimal level 
for L = 8. 

Similar conclusions are verified for the other models, as exemplified in Fig. [2]^a) and 
(b) for the q = 10 Potts model. As in the Ising model, Eq. f lT3|) is not hold for small L's 
(exemplified in part (b) for e = —1.25), which becomes equivalent to (§f)z, for larger (but 
still rather small) L's. It is worth remarking that due to the small number of configurations 
and the discretization of energy, both procedures are not precise in the limit of extreme low 
energies. 

Once established the regime of convergence of methods, we extend the previous analysis 
for the whole range of energy. In Fig. [3] (c), we plot the lni? n as function of n for several 
values of E and L = 10. Note that all curves are linear and cross at (0,0), which gives 

n 

H = for all energies and temperatures, in consistency with results by Beale [14J], where 
the DOS was enumerated for H = 0. In Fig. [3](a) we plot (3l versus the total energy E for 
different L. In order to avoid data overlapping we choose to plot E instead of e = E/V. 
The results for the Ising model show an excellent agreement between estimates of /3l for all 
system sizes (part (6)). In addition, we have also compared (not shown) both schemes at the 
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ferromagnetic-paramagnetic second-order phase transition. The pseudo-critical temperature 
PI may be estimated by the peak in the specific heat C (obtained from the energy numerical 
differentiation). The deviation between (3* L and its asymptotic value f3 c = |ln(l + \/2) 
(obtained here) agrees very well with exact estimates by Beale 

Further, we extend the previous analysis to the Blume-Capel model. In similarity with 
the WL model 15|, numerical simulations were performed for fixed E = Yl(ij) a i°'j an d 
E 2 = In practice, E 2 fixed implies that the number of spins is conserved. The 

one-flip part is restricted to only spins ±1. In Fig. HJa) we plot 0l versus E for different 
E 2 , whereas in the graph (b), we analyze the dependence of /3l on E 2 for E fixed. The 
dependence on L is also showed in Fig. S(c). Note again a very good agreement between 
both approaches, even for small system sizes, supporting once more the equivalence between 
Eqs. ([I]) and (HO]). 

Now we take the Potts model in kind. As in the Ising and BC models, it also presents 
ferromagnetic-paramagnetic phase transitions, exactly located at /3 = ln(l + y/q). For 
q < 4, it is second-order which becomes first-order for q > 4. The case q = 4 presents 
remarkable features, including logarithmic scaling corrections jl], [l7, [3] and a double peak 



probability distribution 



19j, hence an interesting case to be considered. In Figs. Efa) and 



(b) we evaluated /3l for different energies and system sizes (relative small L's but sufficient 
large to imply the validity of Eq. (|13p ). As in the previous examples, we have also found 
an excellent agreement between estimates obtained from microcanonical procedures. In the 
inset of Fig. E£a) we plot the pseudo-critical temperature obtained from the peak in 



the specific 

= (lnL)3/4 

y — L 3/2 



leat C. For q — 4, the deviation of f3* L from its asymptotic value f3 c decays as 



17 



where we found (by using this scaling law) the estimate (5 C = 1.0982(4), 
in excellent agreement with the exact value f3 c = ln(l + = 1.0986123.... 

In the last analysis, we evaluated the microcanonical temperature for the q = 10 and 
q = 20 states Potts model. These are also very interesting cases because, in contrast to all 
previous ones, they possesses discontinuous transitions characterized by a S-like structure, 
a "loop" 16[, hence ideal examples for illustrating the correctness of the present approach. 
Loops for finite systems in the microcanonical ensemble are due to interfacial effects, in 
which the surface tension behaves as 1/L 13|, |l_6|, |2lJ. In contrast, systems simulated in 
the grand-canonical ensemble do not present loops. In Fig. |6] we show the validity of Eq. 
(ITT]) by plotting f} L versus e for q = 10 (part (a)) and q = 20 (part (6)) for L = 20. 
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As in all previous cases, we also have a good agreement between the temperatures, for 
both values of q. However, in contrast with our results, estimates obtained from the WL 
method presents large fluctuations, even using the adaptive windows improvement, taking 
the flatness criterion of 90% for the convergence of H(E) and evaluating the mean DOS over 
10 different seeds. On the other hand, our estimates become more precise as L increases. 
For smaller system sizes they are less accurate (already taking the L for which methods 
are equivalent), despite the accordance with estimates from the WL method. This can be 
understood in the following: Since the number of configurations for fixed E is very large and 
the right side of Eq. (fTTj) is evaluated from only two possible spins, averages become less 
precise for small L. By increasing L, the number of sites with above chosen spins are larger 
and therefore, the averages becomes more precise. This is an interesting point, since on the 
contrary to the WL, under the present approach /3l becomes more precise by increasing L. 
In the inset of each figure, we plot the dependence of the minimum /3£ on L~ l , where the 
straight line have linear coefficients (3q = 1.4267(7) and /?o = 1.702(2), which agrees very 
well with the exact value O = ln(l + V^O) = 1.42606... and /3 = ln(l + V?tt) = 1.699669.... 



IV. CONCLUSIONS 

In this paper we have clarified the fundamental issue of a method proposed at Ref. 
which uses a grand-canonical relationship for calculating microcanonical quantities. The 
study was carried out from the direct comparison with the standard definition of micro- 
canonical temperature. A detailed analysis for three different lattice models yielding first 
and second-order phase transitions sizes was undertaken. Our results show that although 
methods are not equivalent for extreme small system sizes, they converge for relative small 
L's (in practice, our definition becomes equivalent to (|§)i for L = 8). Not only the es- 
timates for finite systems were found to be equivalent, but also the thermodynamic limit 
transition points. The further contribution exploited its advantages. Besides the general- 
ity and easy implementation, thermodynamic quantities are precisely evaluated from rather 
short simulations. Other advantages of the method concerns its low computational cost 
and not requiring a criterion for achieving convergence of results. Once the equivalence has 
been verified in distinct situations, we believe that the present approach may offer a rather 
cheap method for simulating more complex systems, such as lattice models with continuous 
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variables 22[, spin-glasses and polymer systems This will be the subject of ongoing 
work. 
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FIG. 1. Inverse of temperature (3l versus L for different e's for the Ising model. Symbols and 
continuous lines denote the estimates obtained from the present approach and the exact DOS, 
respectively. In {b) we have a mono-log plot of the quantity R n versus n for different L and 
e = —0.75. Curves have been shifted in order to avoid overlapping and continuous lines have been 
used for better visualization of slopes. 
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FIG. 2. Inverse of temperature (3l versus L for different e's for the q = 10 Potts model. Symbols 
and continuous lines denote the estimates obtained from the present approach and the exact DOS, 
respectively. In (6) we have a mono-log plot of the quantity R n versus n for different L and 
e = —1.25. Curves have been shifted in order to avoid overlapping and continuous lines have been 
used for better visualization of slopes. 
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FIG. 3. Inverse of temperature (3l versus total energy E for the Ising model and different system 
sizes L. Symbols and continuous lines (for (a) and (6)) denote the estimates obtained from the 
present approach and the exact DOS, respectively. Inset shows the convergence of methods already 
achieved for L = 8. In (6) we have the (3l versus L for some values of e = E/V . In (c) we have a 
mono-log plot of the quantity R n versus n for L = 10 and different values of E. Slopes of straight 
lines give f3i- In this case, continuous lines have been used for better visualization of slopes. 
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FIG. 4. In (a) and (b) we plot the inverse of temperature (3l as a function of -E (£2) for the BC 
model with E2 (E) fixed and L = 8. Symbols and continuous lines denote estimates obtained from 
the present approach and from Wang-Landau method, respectively. In (c) we show the dependence 
of /3l on L for different values of e and e2 = 2/3. In (d) we plot lnR n versus n for E for e2 = 2/3 
and L = 10. Slopes of straight lines give f3i- In this case, continuous lines have been used for 
better visualization of slopes. 
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FIG. 5. In (a), we plot the inverse of temperature Pl versus E for several L in the q = 4 Potts 
model. Symbols and continuous lines denote estimates obtained from present approach and WL, 
respectively. In (6) we plot the dependence of Pl on L for different values of e. In (c) we plot 
the lni? n versus n for different energies and L = 10. Slopes of straight lines give Pl- In this case, 
continuous lines have been considered for better visualization of slopes. In the inset we plot the 
inverse of the temperature (3 L in which specific heat presents a maximum versus y = ^ n L ^]J ■ The 
straight line has linear coefficient f3 c = 1.0982(4). 
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FIG. 6. Inverse of temperature (3l versus e for the q = 10 Potts model for L = 16 (top) and L = 20 
(bottom). Symbols and continuous lines correspond to the estimates obtained from the present 
approach and WL, respectively. In the inset, we plot the dependence of the minimum f3 L versus 
L , respectively. The straight line has linear coefficient /3q = 1.4267(7). 
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